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Abstract 

The Interstellar Boundary Explorer (IBEX) spacecraft is currently in a highly elliptical orbit around Earth 
with a period near 3:1 resonance with the Moon. Its orbit is oriented so that apogee does not approach the 
Moon. Simulations show this orbit to be remarkably stable over the next twenty years. This article 
examines the dynamics of such orbits in the Circular Restricted 3-Body Problem (CR3BP). We look at 
three types of periodic orbits, each exhibiting a type of symmetry of the CR3BP. For each of the orbit 
types, we assess the local stability using Floquet analysis. Although not all of the periodic solutions are 
stable in the mathematical sense, any divergence is so slow as to produce practical stability over several 
decades. We use Poincare maps with twenty-year propagations to assess the nonlinear stability of the 
orbits, where the perturbation magnitudes are related to the orbit uncertainty for the IBEX mission. 
Finally we show that these orbits belong to a family of orbits connected in a bifurcation diagram that 
exhibits exchange of stability. The analysis of these families of period orbits provides a valuable starting 
point for a mission orbit trade study. 


Introduction 

The Interstellar Boundary Explorer (IBEX) spacecraft is currently in a highly elliptical orbit around Earth 
near 3:1 resonance with the Moon, with spacecraft apogee oriented to stay away from the Moon. As 
described in [1-3], this type of orbit would be useful for space weather missions and shows remarkable 
stability over an interval of several decades. 

This article is an extension of [1], which reviews the trajectory trade study for the IBEX extended 
mission. The goal of this article is to obtain a better understanding of the dynamics of such an orbit, 
especially its stability. Toward this goal we study the dynamics of near 3:1 resonant periodic orbits, in the 
CR3BP for the Earth-Moon system. Poincare recognized the central role of periodic solutions in the 
structure of a dynamical system [2]. As shown in [1], the actual IBEX orbit is quasi-periodic, not 
periodic. However it is common for a family of quasi-periodic orbits to exist near a periodic orbit in a 
Hamiltonian system [3]. 
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There have been several fundamental studies of the dynamics of periodic orbits in the CR3BP, including 
[6-10]. Families of Libration Point Orbits (LPOs) are mapped out in [11-12], and those articles provide 
extensive literature surveys. Resonant periodic orbits have been studied in [4] and [5]. The significance of 
resonant orbits in astrodynamics has long been recognized. The works [14-15] discuss resonances in 
connection with solar system dynamics. Many studies have examined resonance in the dynamics of the 
asteroid belt, including [16-18]. The search for extrasolar planets has identified systems with pairs of 
planets that are near resonance, such as [19-20]. During the Apollo era, Broucke [6] computed numerous 
periodic orbits in the CR3BP of the Earth-Moon system and assessed their stability. There have also been 
extensive studies of cycler orbits that repeatedly approach the two primary bodies in a system such as the 
Earth and Moon [22-23]. However, cycler orbits tend to be unstable due to close approaches with the 
primaries, and so require maneuvers to maintain the periodicity. By contrast, the IBEX orbit requires no 
orbit maintenance maneuvers for at least ten years. The body of work on resonance orbits like IBEX, 
which keeps the spacecraft away from the primary bodies, appears to be small but growing. We did not 
find any orbits in Broucke ’s catalog [6] similar to the ones studies here. Mathews, McGiffin et al. [24-25] 
investigated the use of 2:1 resonance orbits with a lunar gravity assist. These two references were used as 
a basis for the trajectory design of the recently approved Transiting Exoplanet Survey Satellite, to be 
launched in 2017 [7]. Currently Vaquero and Howell [27-29] are investigating the dynamics of resonant 
orbits in the Earth-Moon system. 

The remainder of this paper is structured as follows. 

We first review the orbit properties of the IBEX extended mission. We then summarize the Circular 
Restricted 3-Body Problem (CR3BP) model. In order to compute periodic solutions of the CR3BP, we 
exploit the Mirror Theorem [8]. This approach has been used in previous studies, such as [9], to compute 
LPOs. Based on symmetries in the CR3BP we identify three types of periodic solutions: Planar Mirror, 
Reflection and Axial. To explore the dynamics of near 3:1 resonant solutions, for each of the orbit types 
we compute a particular periodic solution with properties similar to the IBEX orbit. These solutions in the 
CR3BP rotating frame have periods near, but not equal to, the orbit period of the Moon. Because the 
solution period does not match the Moon’s orbit period, the line of apsides of the periodic orbit rotates 
with a secular rate, which can be useful in science mission observations. For each representative solution, 
we use Floquet theory to perform a linear stability analysis. We also show how Lyapunov exponents can 
be computed from Floquet multipliers for a periodic solution. This linear analysis shows that, while the 
solutions may not be stable in the strict mathematical sense, any growth in perturbations over time is so 
slow as to make the solution “long-term stable” [10] compared to a typical spacecraft mission lifetime. 
The Lidov-Kozai mechanism appears to describe a long-term oscillation of orbit inclination and 
eccentricity observed in the quasi-periodic solutions. 

Next we extend the stability analysis to include nonlinear dynamics, using a Poincare section. A similar 
analysis for a Mirror solution was performed by Vaquero and Howell in [27-28]. This numerical 
simulation shows the twenty-year evolution of perturbed solutions, where the perturbation sizes are based 
on the known uncertainties in orbit determination for IBEX. 

While the analysis of three solutions sheds light on the dynamics of near-resonant orbits, the three cases 
we examine may not be entirely representative. To address this concern, we exploit the Cylinder Theorem 
for Hamiltonian systems [11] and compute continuous families of Reflection, Axial and Planar Mirror 
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periodic solutions, and look at the linear stability of each family. Moreover, we show that both the 
Reflection and the Axial branches connect to the Planar Mirror family in a bifurcation diagram. The 
computation of the three branches of periodic solutions can provide a rich starting point for a mission 
trajectory trade study. 

To conclude the analysis we take an initial state from the actual IBEX trajectory, and propagate it for 
twenty years using the CR3BP dynamics. We compare the evolution of the orbit elements under CR3BP 
dynamics with the evolution under a high-fidelity force model to determine what features of the full 
dynamics are captured in the CR3BP dynamics. We conclude with some suggestions for future work. 


IBEX Extended Mission Orbit Dynamics 

The IBEX spacecraft was launched in 2008 into a highly elliptical orbit around the Earth. That orbit 
experienced significant quasiperiodic oscillations in the perigee radius, due to variations in where the 
spacecraft encountered the Moon each orbit. IBEX completed its primary mission in January 2011. The 
project was approved for an extended mission, at which point the project management and science team 
leads decided to move the spacecraft into a more stable orbit. The Keplerian orbit period varies between 
about 8.5 and 9.5 days, with an average of 9.13 days. See Figure 1. The period 9.13 days is very close to 
3:1 resonance with the Earth-Moon sidereal orbit period of 27.3 days. 


SatellHe-Pmdkaon ■ 01 Jun 2013 17:59:55 





Figure 1. Keplerian orbit period evolution for the IBEX orbit, predicted from Jun 2103 for twenty years of the extended 
mission orbit. 

The IBEX orbit, as of April 2013, is inclined to the lunar orbit plane by about 17 degrees. The orbit has an 
apogee radius of about 48 Earth radii (Re) and perigee radius of about 10.5 Re. The Keplerian orbit 
elements vary due to perturbations from lunar and solar gravity, solar radiation pressure and the Earth’s 
nonspherical shape. Because the Moon has a significant influence on the orbit, it is useful to view the 
orbit in the Earth-Moon rotating frame. Due to the near resonance, the orbit appears as a three-leaf pattern 
with apogee near 180 deg or 60 deg from the Moon, as seen in Figure 2. 

The IBEX orbit shown in Figure 2 exhibits a quasiperiodic “nodding” motion in the Earth-Moon rotating 
frame. Figure 1 shows a short-term oscillation corresponding to the orbit period of about 9 days, and a 
second longer-term oscillation with a period of about 9 months. 
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Figure 2. The IBEX orbit viewed in the Earth-Moon rotating frame. In this plot, the Moon is always at the top of the 
circle. The green curve shows the primary mission orbit, where the spacecraft-Earth-Moon angle varied widely at apogee. 
The blue curve shows the extended mission orbit, where the spacecraft-Earth-Moon angle remains near 180 deg or 60 deg 
at apogee. The red curve is the transfer orbit between the primary and extended mission orbits. 


Satellite-Prediction: J2000 Classical Orbit Elements - 02 Jun 2013 18:55:01 
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Figure 3. Twenty -year projection of IBEX orbit Keplerian semimajor axis. A critical feature of the stability of this orbit is 
that the semimajor axis, and so the orbit period, remains steady over a long term, so that the phasing of Moon encounters 
remains fairly constant. 

Figure 3 shows a twenty-year projection of the Keplerian semimajor axis, and exhibits similar oscillations 
to those observed in Figure 1. A key feature of the long-term stability of the IBEX orbit is that the period 
and semimajor axis remain fairly steady, oscillating around a stationary average, which means that the 
orbit apogee maintains regular phasing with respect to the Moon. Figure 4 and Figure 5 show the 
evolution of perigee radius and inclination to the lunar orbit plane. These two plots exhibit the same 10- 
month oscillation together with a much longer-term oscillation with a period near 15 years. 
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Figure 4. Twenty -year projection of IBEX orbit perigee radius, in Earth radii (Re). This curve exhibits two quasiperiodic 
oscillations, one with a period of about 9 months, and another with a period of about 15 years. 


Satellite-Prediction - 01 Jun 2013 17:57:14 



Figure 5. Twenty -year projection of IBEX inclination to the lunar orbit plane, in degrees. This plot exhibits long-term 
trends similar to those of the perigee radius in Figure 4, with a long-term oscillation period of about 15 years. In 
particular, both curves exhibit a minimum about June 2026. 

In April 2013, the argument of perigee was near 270 deg, so the line of apsides was nearly orthogonal to 
the line of nodes. The inertial direction of perigee changed at an average rate of about 1 deg per lunar 
cycle. As noted in [1], this feature is useful in the extended mission because it causes the apogee to move 
in a direction that fills a gap in IBEX’s observations during its initial two-year mission. 

Monte Carlo simulations show that the IBEX extended mission orbit appears stable for decades under 
perturbations on the order of orbit determination uncertainty [10]. This stability is in sharp contrast to the 
primary mission orbit, where perigee radius could not be predicted more than a few years into the future 
[ 1 ]. 


In the remainder of this paper, we study the dynamics of resonant orbits like IBEX in the CR3BP to gain 
an understanding of the orbit dynamics exhibited in Figure 1 to Figure 5. As we shall show, the CR3BP 
model does not completely describe the IBEX orbit dynamics but it does provide valuable insights. 


Page 5 


Circular Restricted 3-Body Problem 

Assume two primary bodies such as the Earth and Moon, each a point mass, orbit each other in circles 
according to Kepler’s laws. A third small body, the spacecraft, is affected gravitationally by the primary 
bodies but does not affect the orbits of the primaries. Mass is scaled so that the total mass of the primaries 
is 1, with the mass of the smaller primary equal to fi. For the Earth-Moon system, /r = 0.0121505. 
Distance is scaled so that the constant distance between the primaries is 1. We assume a distance of 
384,400 km between the Earth and Moon. 

The CR3BP coordinate system rotates with the primaries, and since distance between them is constant, 
the two primaries appear stationary. The origin lies at center of mass of the primaries, and the two 
primaries lie on the x-axis with the larger primary at — /r, the smaller primary at 1 — fi. The z-axis points 
along the orbit normal of the primaries, and the y-axis completes the right-handed orthogonal triad. Time 
is scaled so that the orbit period of the primaries is 2n. For the Earth-Moon system, the CR3BP time 
interval 2n corresponds to 27.28 days. 

The equations of motion of the CR3BP are 

x = v x 

y — Vy 

(1) z — v z 

v' x — 2 v y + x — (1 — ju)(x + n)d^ 3 — nix — 1 + n)d 2 3 
v y = -2v x + y - (1 - \i)ydl 3 - nyd^ 3 
v z — —(1 — /d)zd i 3 — juzdj 3 
where 

d ± = 7 (* + £0 2 + y 2 + z 2 , 
d 2 = 4{x- TT7o^i-y^i-^ 

denote the distances to the primaries. (See [12] for a derivation.) Here prime (‘) denotes a derivative with 
respect to the CR3BP time variable t. The position vector is r = ( x , y, z) T , where superscript T denotes 
the transpose, the velocity vector is v — ( v x , v y , v z ) T , and the state vector for the system is s = 
(x,y,z,v x ,v y ,v z ) T . 


The CR3BP system is autonomous: the equations of motion (1) do not explicitly depend on time. The 
system (1) has one integral of motion, namely the Jacobi integral 

(2) C(s) = 2 U(x, y, z) - (v% + v % + v%) 

where 
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(3) U(x,y,z) = ( x 2 + y 2 )/ 2 + (1 - ii)^ 1 + ^d 2 1 ■ 


The conservation of the Jacobi integral is associated with the Tisserand criterion, which asserts that 
quantity 

(4) J a (i - e 2 )cos(j) 

is nearly constant. (See, for example, [13] , p. 72.) In expression (4), a is the semimajor axis of the orbit, 
e is the eccentricity, and i is the inclination to the orbit plane of the primaries. We note in the Appendix 
that the system (1) is a Hamiltonian system, though not in standard form, where the Hamiltonian H = 
-C/2. 

In this paper we focus on the computation and analysis of periodic orbits in the CR3BP. The CR3BP has 
symmetries that can be exploited in the computation of periodic orbits, and a fundamental tool is the 
Mirror Theorem in [8]. A mirror configuration in the N-body problem is a configuration where each 
velocity is orthogonal to each position. For the CR3BP there are three important mirror configurations: 

• Planar Mirror: The particle lies at position (x, 0,0) on the x-axis, with velocity (0, v y , 0) 
orthogonal to the x-axis and in the x — y plane. 

• Reflection: The particle lies at position (x, 0, z) on the x — z plane with velocity (0, v yt 0) 
orthogonal to the x — z plane. 

• Axial: The particle lies at position (x, 0,0) on the x-axis with velocity (0, v y , v z ) orthogonal to 
the x-axis. 

The three resonant families are analogous to more familiar families of LPOs that exhibit the same 
symmetries [14] [15]: 

• Planar Mirror resonant orbits are analogous to Lyapunov LPOs 

• Mirror resonant orbits are analogous to Halo LPOs 

• Axial resonant orbits are analogous to Axial LPOs 

It is shown in the section on Families of Near 3:1 Resonant Orbits that these analogies between resonant 
orbits and LPOs is also exhibited in the connections between the families in a bifurcation diagram. 

The Planar Mirror configuration is actually a special case of both the Reflection and the Axial mirror 
configurations. The Mirror Theorem asserts that if the solution is in a mirror configuration at time t — 0, 
and in another mirror configuration at t — T /2, with T > 0, then the solution is periodic with period T. 
The Mirror Theorem has been exploited in works such as [9] to compute periodic orbits in the CR3BP. 
For each of the mirror configurations above we compute a type of periodic orbit: 

• A Planar Mirror periodic orbit has Planar mirror configurations at its initial point and midpoint. A 
Planar Mirror periodic orbit is symmetric both under reflection across the x — z plane and under 
rotation by 180 degrees about the x-axis. 
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• A Reflection periodic orbit has Reflection mirror configurations at its initial point and midpoint. 
A Reflection periodic orbit is so named because it is symmetric under reflection across the x — z 
plane. 

• An Axial periodic orbit has Axial mirror configurations at its initial point and midpoint. An Axial 
periodic orbit is so named because it is symmetric under rotation by 180 degrees about the x-axis. 

In this paper the solutions are computed in the rotating CR3BP frame. However for spacecraft mission 
analysis it is essential to understand the solution properties in an inertial frame: its orientation, its shape, 
and its evolution. For this purpose we have defined an inertial “EM” frame (for Earth-Moon) with the 
following properties: The origin of the EM frame is the Earth. The Moon orbits the Earth in a circle in the 
EM plane. At time t = 0, the axes of the inertial EM frame and the rotating CR3BP frame are aligned. 
Thus at t = 0, the Moon lies along the +x axis of the inertial EM frame. 


Three Representative Orbits 

No single periodic orbit in the CR3BP fully represents the actual IBEX orbit. To gain some understanding 
of the dynamics of this class of orbits, we first consider a representative example of each of the three 
types of periodic solutions described above, and then later look at the families of periodic orbits. To 
compute each periodic solution we employ a differential correction shooting method based on the Mirror 
Theorem. The shooting method process is adapted to the symmetries of each case. For a Planar Mirror 
case, the shooting method fixes the initial values of y, z, v x and v z at 0, and also fixes the initial value of 
x. We choose the initial value of x to achieve a desired apogee radius. The estimated period 7 and the 
initial value of v y are allowed to vary to achieve the goals that at time t = 7/ 2, the values of y and v x 
equal 0. After each step of the differential correction, the estimate of the period 7 is improved for the 
updated state by determining the time 7 /2 when the trajectory crosses the y = 0 plane. 

For the Reflection case, the shooting method is similar to that for the Planar Mirror case but must also 
handle the z component. The method fixes the initial values of y, v x and v z at 0, and also fixes the initial 
orbit inclination i 0 to the EM plane. The estimated period 7 together with the initial values of x and v y 
are allowed to vary to achieve the goals that at time t = 7/ 2, the values of y, v x and v z equal 0. The 
initial value of z is chosen to achieve the initial inclination, given the other initial state components. After 
each step of the differential correction, the estimate of the period 7 is improved for the updated state by 
determining the time 7/2 when the trajectory crosses the y = 0 plane. 

For the Axial case, the shooting method fixes the initial values of y, z and v x at 0, and fixes the initial 
orbit inclination i 0 . The estimated period 7 together with the initial values of x and v y are allowed to vary 
in order to achieve the goals that at time t = 7/2, the values of y, z and v x equal 0. The initial value of v z 
is chosen to achieve the initial inclination, given the other initial state components. After each step of the 
differential correction, the estimate of the period 7 is improved for the updated state by determining the 
time 7 /2 when the trajectory crosses the v x = 0 hyperplane. We focus on hyperplane v x = 0 for the 
Axial case based on its symmetry. 


Planar Mirror Case 

Figure 6 shows two plots of a Planar Mirror solution for the case where the apogee radius was fixed at 
49.95 Earth radii (Re), similar to the IBEX orbit when we performed this analysis in January 2012. The 
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left plot shows the view in the CR3BP frame while the right plot shows the view in the EM inertial frame. 
The reflection symmetry of the orbit across the CR3BP x — z plane is apparent in the left plot. The 
components of the initial state, together with some key parameters, are shown in Table 1. The plot in the 
CR3BP frame shares the three-leaf structure of the actual IBEX orbit in Figure 2. 


CR3BP Periodic Solution in CR3BP Frame CR3BP Periodic Solution in Inertial EM Frame 




Figure 6. Views of the Planar Mirror periodic orbit. On the left in the solution in the CR3BP frame. Note the mirror 
symmetry across the x axis. On the right is the orbit in the inertial EM frame. Because the orbit period is not the same as 
the orbit period of the primary bodies, the orbit in the inertial frame is not periodic. The inertial plot shows evolution 
over one year. 

During the design of the IBEX orbit, one goal was to have the line of apsides rotate in a desired direction, 
so that the observations would cover a direction in space not covered during the initial two-year mission. 
It is apparent in the right plot in Figure 6 that the line of apsides rotates in inertial space. A critical factor 
in the apsidal rotation is the difference between the orbit period T in the CR3BP time unit and the orbit 
period 2n of the primary bodies. Figure 7 illustrates this concept. For example, the Planar Mirror orbit 
considered here has an orbit period of T = 6.2746, slightly less than 2n. After one satellite orbit period T 
in the CR3BP frame, the satellite is again at apogee and is in opposition with the Moon. Now consider 
what happens after one orbit period T in the inertial frame. In Figure 7, the spacecraft and Moon are both 
located on the x-axis. After one orbit period T , the spacecraft goes from one apogee to the next, and the 
spacecraft and Moon are again in opposition. Poincare makes a similar observation in [2], Vol. 1, section 
39. Because the period T is less than 2n , the Moon does not complete an entire orbit. As a result, the line 
of apsides rotates clockwise. The average apsidal rate is T — 2n radians per lunar cycle. To measure 
apsidal direction, we define the perigee angle to be the angle from the x-axis to the perigee direction 
projected into the EM x — y plane. Figure 8 shows the evolution of the perigee angle over time for the 
Planar Mirror solution computed here. 
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Figure 7. Rotation of the line of apsides in the EM inertial frame for a nearly 3:1 resonant orbit. The Moon (marked by a 
small circle) follows the black outer orbit, and the spacecraft (marked by an asterisk) follows the blue orbit. At time t = 0 
the spacecraft is at apogee and is in opposition with the Moon, with both on the x-axis. After one orbit period T in the 
CR3BP frame, or three orbits in inertial space, the spacecraft is once again at apogee and in opposition with the Moon. In 
this illustration we assume the period T is less than the Moon’s orbit period. Therefore after time T the Moon has not 
completed a full orbit, so does not lie on the x-axis. Because the spacecraft is in opposition with the Moon at time T, the 
spacecraft line of apsides must have rotated from its original direction along the x-axis. The angle of rotation per orbit is 
proportional to the difference between the Moon orbit period and the spacecraft orbit period T. (The angle of rotation is 
exaggerated for illustration. It would typically be at most a few degrees per period T.) 


Perigee Angle vs. time 



Figure 8. Perigee angle versus time for a near 3:1 EM-resonant orbit over one year. This plot is for a planar orbit with 
radius of perigee near 8.15 Re which has a period in the CR3BP frame equal to 6.2746. Therefore after one orbit period 
the Moon is about 0.6 deg short of completing one orbit, and after one year, the perigee angle rotates by about 6.6 deg. In 
addition to the secular variation, the perigee angle also oscillates each orbit period as the spacecraft-Moon geometry 
varies. 
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Reflection Case 

Figure 9 and Figure 10 shows a Reflection periodic solution with initial inclination of 23 deg to the EM 
plane, which is about the average inclination shown in Figure 5 for the predicted IBEX orbit. The initial 
state components and other parameters are again in Table 1. The CR3BP plot of this solution in Figure 9 
exhibits the same three-leaf shape as in Figure 6. However the perigee radius for the Reflection solution is 
more than 10 R E , higher than the perigee radius for the Planar Mirror orbit. Our computations indicate 
that there is one Reflection solution for each initial inclination. Figure 9 shows the projection of the orbit 
into the x — y plane, while Figure 10 shows a skew view in three dimensions. In the CR3BP frame, the 
near 3:1 resonant orbit looks similar to a three-petal flower. In the inertial EM frame, the initial line of 
nodes is along the y-axis. 



Figure 9. Views of the Reflection periodic orbit in CR3BP and inertial EM frame, projected into the xy-plane. 


CR3BP Periodic Solution in CR3BP Frame 


CR3BP Periodic Solution in Inertial EM Frame 



Figure 10. Skew views of the Reflection periodic orbit in CR3BP and inertial EM frame, projected into the xy-plane. The 
inertial EM frame on the right emphasizes the inclination of the spacecraft orbit to the Moon orbit. 
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Axial Case 

Figure 1 1 and Figure 12 show the Axial solution for initial inclination 23 deg. As shown in the skew view 
in Figure 12, an Axial solution in the CR3BP frame has a shape like a three-blade propeller. It can be 
difficult to see in a single skew view, but the Axial solution is symmetric under a rotation by 180 deg 
about the x-axis. (The 3:1 resonant solution is also symmetric under a rotation by 120 deg about the z- 
axis. However the x-axis symmetry is a fundamental symmetry of the dynamics.) The view in the EM 
inertial frame shows that the initial line of nodes lies along the x-axis. 


CR3BP Periodic Solution in CR3BP Frame CR3BP Periodic Solution in Inertial EM Frame 




Figure 11. Views of the Axial solution, projected into xy plane 



Figure 12. Skew views of the Axial solution. In the CR3BP frame on left, the 3:1 resonant Axial solution resembles a 3- 
blade propeller. In the EM frame on right, the orbit is inclined to the xy plane with the apsides essentially along the x- 
axis. 


Local Stability Analysis 

This study was undertaken to comprehend the apparent stability of the near-resonant orbits like the IBEX 
orbit. There are numerous types of stability defined in dynamical systems, including asymptotic, linear, 
orbital and spectral stability. See, for example, [11] and [34-35]. To study the local stability of a solution 
we linearize the dynamics around the solution. The CR3BP is autonomous, and the solution is periodic 
with period T , so the linearization about the reference periodic orbit is periodic with period T. Floquet 
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theory shows that, to assess the stability of the reference orbit, it is sufficient to examine the monodromy 
matrix M, which is the state transition matrix 0(7) for the time interval T [16]. 



Planar Mirror Solution 

Reflection Solution 

Axial Solution 

CR3BP x 

-0.84094370 

-0.74257943 

-0.85249293 

CR3BP y 

0 

0 

0 

CR3BP z 

0 

0.31004867 

0 

CR3BP v x 

0 

0 

0 

CR3BP v y 

0.25043310 

0.06950217 

0.33447503 

CR3BP v z 

0 

0 

-0.21472796 

CR3BP period 

6.27459754 

6.26105536 

6.27964378 

Period (days) 

27.24731418 

27.18850753 

27.26922740 

Jacobi integral 

3.04158388 

3.05024235 

2.93303141 

Initial Apogee Radius (Re) 

49.95 

47.82 

50.64 

Initial Perigee Radius (deg) 

8.15 

10.18 

7.46 

Initial inclination (deg) 

0.00 

23.00 

23.00 

Floquet Multiplier 1 

1.0000 

1.0000 

1.0000 

Floquet Multiplier 2 

1.0000 

1.0000 

1.0000 

Floquet Multiplier 3 

0.7172 - 0.6968i 

0.8097 - 0.5868i 

0.7194 - 0.6946i 

Floquet Multiplier 4 

0.7172 + 0.6968i 

0.8097 + 0.5868i 

0.7194 + 0.6946i 

Floquet Multiplier 5 

0.9922 

0.9997 - 0.0227i 

0.9907 

Floquet Multiplier 6 

1.0078 

0.9997 + 0.0227i 

1.0094 

Period of first 
quasiperiodic mode Q1 

(days) 

222.0 

272.4 

223.1 

Period of second 
quasiperiodic mode Q2 (days) 


6167 


Time for unstable mode U to 
grow by factor of 10 (days) 

8090 


6719 


Table 1. Initial states for three representative periodic orbits, with additional dynamical data. Each solution has a pair of 
Floquet multipliers equal to +1, and a complex conjugate pair with magnitude 1 corresponding to a quasiperiodic mode. 
The Reflection solution has a second complex conjugate pair, while the Planar Mirror and the Axial each have a real 
inverse pair. For each quasiperiodic pair, the corresponding period is shown, whereas for each unstable mode, the time to 
grow by a factor of 0 is shown. In all three cases, the final pair is very close to 1, yielding either a very long period 
oscillation or a very slow growing instability. 

Because the CR3BP system is Hamiltonian, the monodromy matrix M has some special properties, 
discussed in the Appendix. For stability analysis the eigenvalues of M, called the Floquet multipliers, play 
a critical role. The enumeration of possible configurations of Floquet multipliers in the CR3BP is 
described in [17] (p. 157 and p. 173), [16] and [18]. In this study we encounter three configurations: 


• A pair of values +1 

• A complex conjugate pair A and A, each with magnitude 1 

• A real pair of inverses A and 1/A, so the product is 1. 

In an autonomous Hamiltonian system the determinant of M, which is the product of eigenvalues, must 
equal 1 . For a periodic solution in a Hamiltonian system, there must be a pair of Floquet multipliers equal 
to +1. Moreover, as noted in the Appendix, one eigenvector corresponding to eigenvalue +1 is the right 
hand side of the differential equation (1) evaluated at the initial state. Floquet multipliers measure how 
perturbation components called Floquet modes evolve. If a Floquet multiplier A corresponds to 
eigenvector 17, and the periodic solution is perturbed by the amount £ v , then linear analysis predicts that 
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after m orbit periods the perturbation becomes A m £ v. Thus if a Floquet multiplier has magnitude greater 
than 1, then that Floquet mode grows exponentially, and the solution is unstable. However, among the 
solutions studied here, when there does exist an unstable Floquet multiplier then it is only slightly larger 
than 1, and the growth in the unstable direction can take decades to reach a significant size. In an 
autonomous Hamiltonian system, the Hamiltonian is conserved and phase space volume is conserved, so 
a solution cannot be asymptotically stable. If all of the Floquet multipliers have magnitude 1, then the 
solution is called spectrally stable. If in addition there is a complete set of eigenvectors, then the solution 
is linearly stable [11]. 

In the deficient case where the pair of eigenvalues +1 have only one eigenvector, the Jordan canonical 
form of M is not diagonalizable and contains a 2-by-2 block with diagonals +1, and there is an associated 
generalized eigenvector. Weisel and Pohlen [18] showed that this situation occurs when the periodic 
solution is part of a continuous family of periodic solutions, and the period varies with an integral of 
motion. A familiar example of deficiency is Keplerian motion where the orbit period depends on the 
semimajor axis. In this situation, the periodic solution itself is not stable. (Suppose that a spacecraft is 
placed in a circular orbit, and a second spacecraft is placed in a slightly higher circular orbit near the first 
satellite. Because the second satellite has a longer orbit period it drifts away from the first satellite in the 
along-track direction, initially at a linear rate.) This is a very mild sort of instability, where the separation 
rate is only linear in time. Although the two solutions drift apart, their two orbits may remain close 
together, a phenomenon called orbital stability. (See [19], p. 274.) 

Figure 27 shows that each of the periodic solutions studied here is in fact part of a family where the 
period varies with Jacobi integral, so the Jordan block with eigenvalue +1 is not diagonalizable. For each 
periodic solution we computed the Jordan canonical form of M using the Matlab function eig. However, 
the numerical errors inherent in finite-precision calculations make it challenging to detect eigenvector 
degeneracy. Instead eig produces a pair of distinct, approximate eigenvalues near +1. The Appendix 
describes a way to use the Matlab results to obtain the Jordan canonical form. 

In the eigenvalue decomposition of each periodic solution we also find at least one complex conjugate 
pair with magnitude 1, corresponding to a quasiperiodic mode. For such a pair it is useful to convert to a 
real Jordan block form in a standard way [20]: A pair of eigenvalues a ± i/3 in the complex Jordan 
canonical form yields a 2-by-2 real Jordan block of the form 


The corresponding vectors, the Floquet modes, are the real and imaginary parts of the eigenvectors. 


See, for example [21]. 

The result of the Floquet analysis for each of the three types of periodic solutions is summarized in Table 
1. The three solutions have similar Floquet structures. 

• A pair of multipliers +1, with modes labeled PI a and Plb: 



Floquet modes can be used to control motion around a periodic solution based on the natural dynamics. 


Page 14 


o PI a lies along the eigenvector which is the time derivative of the state, given by equation 
(1). This corresponds to translation along the periodic solution, so the perturbation does 
not grow over time. 

o Plb lies along the generalized eigenvector. Perturbation along this direction produces a 
small change in Jacobi integral and a slightly different period. The perturbed orbit is 
nearly the same as the nominal periodic solution. However the perturbed solution drifts 
away from the periodic solution linearly in time. 

• A complex conjugate pair of multipliers with magnitude 1, with modes Qla and Qlb: 

o This produces a quasiperiodic oscillation about the periodic solution. A complex 
multiplier A = a + i/3 yields a rotation of angle 6 = atanl^p, a) radians about the 
periodic solution each orbit period. The oscillation period is therefore 2nT/9 . For the 
Planar Mirror, Axial and Reflection solutions, the oscillation periods for this mode are 
222, 223 and 262 days, respectively. 

The Reflection solution has a second pair of quasiperiodic modes Q2a and Q2b with a very long period of 
6167 days or 16.88 years. Thus the Reflection solution is spectrally stable. In contrast, the Planar Mirror 
and the Axial solutions have a stable mode S and an unstable mode U corresponding to a real inverse pair 
of Floquet multipliers close to 1. For the Planar Mirror solution, the stable and unstable modes are both in 
the z direction, so that in-plane motion is spectrally stable. Despite this distinction between the long-term 
dynamics of the Reflection solution versus the Planar and Axial solutions, the Floquet multipliers in the 
third pair for all of these solutions are all close to +1. Numerical results below show that a very long- 
period oscillation and a slowly growing unstable mode appear very similar for many years. 

Floquet Modes of the Planar Mirror Solution 

To understand more deeply the dynamics near each solution, we applied perturbations in each of the 
Floquet modes and propagated each perturbed solution for twenty years. As a reference we used a 
perturbation scale of o — le — 5 in the CR3BP components. This corresponds to perturbations of 3.8 km 
in position, and le-5 km/sec in velocity. These values are close to the standard deviations in orbit 
determination solutions for IBEX at apogee. For example, Figure 13 shows the quasiperiodic solution 
produced when we perturb the Planar solution by 10,000 o along the quasiperiodic Floquet mode Ql. 

This quasiperiodic solution looks similar to the IBEX orbit shown in Figure 2. 
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Quasiperiodic Solution about Planar Solution 



Figure 13. Quasiperiodic solution produced by perturbing the Planar Mirror periodic solution along the first oscillatory 
Q1 mode. This quasiperiodic solution is obtained from a perturbation in the Q1 component with magnitude 10,000 o. 

In most of our simulations, the perturbation magnitude is 1000 <7, except where noted. Figure 14 shows 
the variations in semimajor axis for the various modes. The right plot in Figure 14 clearly exhibits the 
quasiperiodic oscillations with period around 222 days. The left plot shows that in the unstable mode (red 
curve), the oscillation in the semimajor axis grows slightly over twenty years from about 25 km to 50 km. 
However, even in the unstable mode, the average value of semimajor axis remains the same, so that the 
phasing of encounters with the Moon remains about the same. By contrast, Figure 15 and Figure 16 show 
that the unstable mode does cause the perigee radius and inclination to change over time. For the Planar 
solution, Table 1 shows that the largest Floquet multiplier is 1.0078. We propagate for 20 years or 267.6 
lunar cycles. Thus we would expect the perturbation in the unstable mode to multiply by a factor of 
1.0078 267,6 = 8.00, and the stable mode to multiply by a factor of 1/8.00 = 0.125. Figure 15 is 
consistent with those values. There is a much slower growth in the oscillation of the semimajor axis, 
shown in Figure 14. These plots show the evolution of a perturbation in the individual Floquet modes. A 
random perturbation in the initial state would be a linear combination of the Floquet modes. Thus, only a 
portion of the initial perturbation would be in the unstable direction. Moreover, another component of the 
perturbation would be in the stable direction, and its decay would mitigate the growth due to the unstable 
direction. Thus we would expect the typical growth rate to be slower than the unstable mode alone would 
predict, at least until the stable component has decayed to near 0. 
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Planar: SMA Variation vs. Time for Floquet Modes 
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Planar: SMA Variation vs. Time for Quasiperiodic Mode Q1 



Figure 14. Variation in semimajor axis for Planar solution, under perturbations in the Floquet modes with magnitude 
1000 a. The plot on the right shows the evolution of mode Q1 for the first three years, and more clearly exhibits the 
period near 220 days. 


Planar: Perigee Radius Perturbation vs. Time for Floquet Modes Planar: Inclination Perturbation vs. Time for Floquet Modes 




Figure 15. Variations in perigee radius (left) and inclination (right) for Planar solution, for perturbations with magnitude 
1000 <t. The unstable mode causes the most significant change in perigee radius. For inclination, only the stable and 
unstable modes are shown since the other modes do not affect out-of-plane motion. 

To gain a better understanding of the effect of perturbation magnitude, we ran simulations with 
magnitudes 100 a, 1000 a and 10,000 a. Some of the results are shown in Figure 16. For magnitudes 
100 crand 1000 awe observe exponential growth in perigee radius and in inclination for the twenty-year 
duration of the simulation. It is noteworthy that for amplitude 10,000 a, the growth in perigee radius and 
in inclination perturbation appears exponential at first, but then levels off after about 16 years. This 
behavior may be associated with the Tisserand criterion (4), which together with the preservation of 
semimajor axis puts bounds on the range of eccentricity and inclination values. It was also found that if 
we apply a perturbation magnitude of about 22,000 a in the quasiperiodic mode Ql, the oscillations like 
those seen in Figure 13 become so large that the orbit apogee gets close to the Moon, disrupting the 
resonance and producing a chaotic orbit. 
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Planar: Perigee Radius Perturbation vs. Time for Unstable Mode 




Figure 16. Variations in perigee radius (left) and inclination (right) for different perturbation magnitudes. The vertical 
axis in each plot uses a logarithmic scale, so a linear trend in the curve indicates exponential growth. For perturbation 
magnitude 10,000 (7, both perigee radius and inclination appear to level off and remain bounded. 

Figures 14, 15 and 16 have several characteristics similar to the projected IBEX orbit shown in Figures 3, 
4 and 5, with some notable differences. The semimajor axis maintains an essentially constant average in 
both cases. The perigee radius and the inclination show an oscillation on the order of months, although 
the oscillation period shown for IBEX is slightly longer than the 222-day oscillation period for the Planar 
solution. The perturbed Planar solution also shows long-term growth, on the time scale of ten to twenty 
years. For initial perturbations on the scale of 10,000 a, the perturbation remains bounded, just as the 
IBEX oscillations remain bounded in this time frame. 

Floquet Modes of the Reflection Solution 

In this section we examine the Floquet modes for the Reflection solution, using much the same analysis 
as for the Planar solution. The key distinction is that the Reflection solution has a second quasiperiodic 
mode Q2, with a period of 16.88 years, instead of the stable/unstable mode pair for the Planar solution. 
The second quasiperiodic mode is apparent in Figures 17 and 18 plotting the evolution of perigee radius 
and inclination. Even when the perturbation magnitude is increased to 10,000 a in Figure 19, the long- 
term oscillation persists. The evolution of orbit elements for the Reflection solution, with the three 
periods on the scale of 9 days, 10 months and 17 years, is qualitatively very similar to the IBEX orbit 
evolution, even though the two orbits have significantly different initial perigee radius values. 
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Figure 17. Variation in perigee radius for the Reflection solution. The plot on the left shows that the quasiperiodic mode 
Q2 oscillates with period near 17 years, as predicted. The plot on the right focuses on the initial 3-year interval, and 
shows more clearly the variations due to mode Q1 with period near 270 days. 


Reflection: Inclination Perturbation vs. Time for Floquet Modes 



Figure 18. Variations in inclination for Reflection solution. Note the Q2 variations in particular, which appears to be in 
phase with the perigee radius in Figure 17. 


Reflection: Inclination Perturbation vs. Time for Q1 Mode 



Figure 19. Variation in inclination for Reflection solution, for various perturbation magnitudes. As in Figure 16, we use a 
logarithmic scale to capture the wide range of values. Although the initial perturbation is in the Q1 mode, the oscillation 
for magnitude 10,000 a exhibits a 17-year period, like the Q2 mode. 
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Floquet Modes of Axial Solution 

The Floquet modes of the Axial solutions are similar to those for the Planar solution, but there are some 
distinct differences. In Figure 20, the average semimajor axis value stays essential constant, as with the 
Planar case. However, for the Axial unstable mode, the growth in the oscillation magnitude is closer to 
the factor of 1.0094 267,6 = 12.23 predicted by linear theory. 


Axial: SMA Perturbation vs. Time for Floquet Modes 
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Figure 20. Variation is semimajor axis for Axial solution. Note in particular the growth in the oscillation magnitude for 
the unstable mode, although the average value remains nearly constant. 

Figure 21 shows that for the Axial solution, the unstable mode causes both the perigee radius and the 
inclination to decrease, in contrast to the Planar solution where both increase. The two solutions share the 
nonlinear behavior that for an initial perturbation magnitude of 10,000 a, the perturbation initially grows 
then remains bounded. In Figure 22, the perturbation in the U mode for the Axial solution appears to 
behave as a long-period oscillation. That behavior is similar to what we see for the IBEX orbit in Figure 4 
and Figure 5. 


Axial: Perigee Radius Perturbation vs. Time for Floquet Modes 



Axial: Inclination Perturbation vs. Time for Floquet Modes 
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Figure 21. Variation in perigee radius (left) and inclination (right) for Axial solution. 
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Axial: Perigee Radius Perturbation vs. Time for Unstable Mode 



Axial: Inclination Perturbation vs. Time for Unstable Mode 



Figure 22. Variation in perigee radius (left) and inclination (right) for Axial solution, for various perturbation 
magnitudes. For the magnitude 10,000 sigma, the perturbation appears to remain bounded over decades. 

Another way to measure the time rate of change of perturbations, closely related to the Floquet multiplier, 
is the Lyapunov exponent. A Lyapunov exponent is a more general tool than a Floquet multiplier, because 
a Lyapunov exponent can be computed for a nonperiodic solution. Let {ray (t)), j = 1 to 6, be the 
eigenvalues of the state transition matrix s 0 ) for a given initial state s 0 Then the Lyapunov 
exponents Ay are defined by 

1 

(5) h = hm-Zn(ra 7 (t)) 

7 £-» oo t v 7 7 

when the limit exist. (See [22], p. 67.) For a periodic solution with initial state s 0 and period T , let {/iy} be 
the Floquet multipliers, i.e. the eigenvalues of the monodromy matrix M = <P(T, s 0 ). It follows that 
M fe = «D(/cT, s 0 ), so rrij(kT ) = [ij for k = 0,1,2, ... We can use these observations in equation (5) to 
show that the Lyapunov exponents for a periodic solution are related to the Floquet multipliers by 

(6) /jLj = exp[XjT\j = 1 to 6 

(See [22], p. 68.) In this study we used the Lyapunov exponents associated with the quasiperiodic modes 
to compute the oscillation periods. 


Nonlinear Stability Assessment using Poincare Maps 

The stability analysis to this point has been based on the linearized Floquet analysis. For a more thorough 
assessment of the stability we must consider a variety of perturbations. The Poincare map is a tool that, 
like Floquet analysis, allows us to reduce the dynamics to a discrete-time mapping [22]. In a Poincare 
map we first define a Poincare section, which is a five-dimensional hypersurface in the six-dimensional 
state space. For the map we record the discrete points where each trajectory crosses the Poincare section. 
Like the differential correction methods described above, we selected the Poincare section to suit the 
periodic solution. For the Planar Mirror and the Reflection solutions, we used the hyperplane y = 0 as the 
Poincare section. For the Axial solution, we used the hyperplane v x = 0. We then generated a set of 
random perturbations from the initial state of the periodic solution, propagated the perturbed solution for 
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twenty years, and recorded the Poincare section crossings. Each of the six components of each 
perturbation was a random draw from a normal distribution with mean zero and a fixed standard 
deviation: 150 draws with standard deviation lOOcr, and 150 draws with standard deviation lOOOcr. 
Because we perturb each of the six components, the standard deviation of the state vector perturbation 
magnitude is V6 = 2.45 times the standard deviation of each component perturbation. It is common with 
Poincare maps to restrict the perturbed states to the same value of Jacobi integral as the nominal solution, 
and we have done so here. For each random draw we computed the normal to the level set of the Jacobi 
integral, and iteratively projected the perturbed solution toward the level set of the nominal periodic 
solution. This “isosurface projection” only required a few iterations for each draw of the initial 
perturbation. 


Poincare Section for Planar Mirror Solution 



Poincare Section for Planar Mirror Solution 



CR3BP x 


Figure 23. Poincare map for the Planar Mirror solution. The Poincare section for this mapping is the hyperplane y = 0, 
and the resulting map is plotted in the x — v x plane. The black curve represents the Planar Mirror periodic solution. The 
green dots represent the evolution of perturbations with standard deviation 100<r, while the blue dots represent the 
evolution of perturbations with standard 1000<r. The distribution of the data points is due primarily to quasi-periodic 
oscillations. The plot on the right focuses on the region near the center of the left plot. 

Figure 23 shows the computed Poincare map for the Planar Mirror solution, projected into the x — v x 
plane. A similar mapping appears in [23]. The black curve shows the Planar Mirror solution, the green 
dots show the crossings for standard deviation lOOcr and the blue dots show the crossings for standard 
deviation lOOOcr. The results indicate stability under large perturbations: The Poincare map points remain 
near the nominal solution. The spread of the crossing points is due to the quasiperiodic oscillations. The 
plot on the right in Figure 23 shows details around the crossing at the middle of the left plot. The ring-like 
patterns are indicative of the quasiperiodic oscillations. 

Figure 24 shows the Poincare map for the Reflection solution. The patterns are similar to those in Figure 
23, but more complex due to the out-of-plane motion. Figure 25 shows the Poincare map for the Axial 
solution, projected in the x — y plane, and again indicates effective stability under the perturbations. 
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Poincare Section for Reflection Solution 



Poincare Section for Reflection Solution 
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Figure 24. Poincare map for the Reflection solution. 


Poincare Section for Axial Solution 



Poincare Section for Axial Solution 
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Figure 25. Poincare map for the Axial solution. The Poincare section for this mapping is the hyperplane v x = 0, and the 
map is plotted in the x — y plane. The shape of the region shown in the right plot, with a “spike” protruding to the 
downward left, suggests that there are two overlapping patterns corresponding to distinct oscillations. 


Families of Near 3:1 Resonant Orbits 

Thus far we have examined in detail three particular periodic solutions similar to the IBEX orbit. For a 
deeper understanding of the dynamics we looked at other nearby periodic solutions. A fundamental result 
in the study of periodic solutions in a Hamiltonian system is the Cylinder Theorem [11], Theorem 6.5.1, 
p. 82, which states that in generic circumstances a periodic solution lies on a smooth cylinder of periodic 
solutions parameterized by an integral of motion. Using elementary continuation methods, we generated 
portions of the Planar Mirror, Reflection and Axial families starting from the particular solutions studied 
above. Figure 26 shows a projection of these families in a plot of perigee radius vs. inclination. In Figure 
26, each point on the three curves represents the initial state of a periodic solution. The particular Planar 
solution considered above has an initial perigee radius of 8.15 Re, while the particular Reflection and 
Axial solutions are at initial inclination 23 deg. The Reflection family was computed from inclination 0 
up to 67 deg, where the perigee radius went below 2 Re. For the Axial family we computed the branch 
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from 0 deg to 81.1 deg inclination, at which point the orbit becomes nearly circular and the initial velocity 
in the rotating frame becomes nearly vertical. Although the velocity in the rotating frame is nearly 
vertical, the velocity in the inertial frame is not because the transformation from rotating to inertial frame 
introduces the term ft) xr in the velocity. This explains why the inclination in the inertial frame is 81.1 
deg, not 90 deg. Figure 26 represents three key families of periodic orbits with apogee away from the 
Moon, but it is not exhaustive. In particular, the family of Reflection orbits we computed is the branch 
with perigee above the x — y plane, analogous to the “Northern” branch of halo orbit [9]. There is a 
symmetric “Southern” branch with perigee below the x — y plane. Similarly the branch of Axial orbits we 
computed has perigee at the descending node, but there is a symmetric branch with perigee at the 
ascending node. 


3:1 Periodic Families: Perigee Radius vs. Inclination 



Figure 26. A bifurcation diagram of the families of periodic solutions. A point on each curve represents the initial point of 
a periodic orbit. The vertical blue line at inclination = 0 represents the Planar orbits with various Rp values. The 
Reflection orbit family and Axial orbit family each intersect the Planar orbit family, at Rp = 11.64 Re and at Rp = 6.85 
Re, respectively. Although the Reflection and Axial families appear to intersect each other in this two-dimensional 
projection, they do not actually meet. The S and U labels indicate whether the solutions on a given segment are stable or 
unstable in the spectral sense. 

A significant feature of this diagram is that the Reflection and Axial families each bifurcate from the 
Planar Mirror family, at different points. The initial states and dynamic properties of the bifurcation 
points are given in Table 2. Note that each bifurcation point in Table 2 has two pairs of Floquet 
multipliers +1, showing that it belongs to two families of periodic solutions. In Figure 26 we have labeled 
segments of the branches with an S or U to indicate respectively whether it is stable or unstable in the 
spectral sense. In particular, the Planar Mirror solutions are spectrally unstable for perigee radius between 
6.85 and 1 1.64 Re, but spectrally stable outside that range. It is not uncommon for an exchange of 
stability to occur at bifurcation points [24], as pairs of eigenvalues merge at +1 then split. 

Figure 27 provides another view of the bifurcation diagram, showing the CR3BP orbit period versus 
Jacobi integral. This diagram shows that the orbit period does vary with the Jacobi integral for each 
family, so the Floquet multiplier +1 must occur in a 2-by-2 block in the Jordan canonical form [18]. 
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3:1 Periodic Families: Period vs. Jacobi Integral 



Figure 27. Another view of the bifurcation diagram for periodic families. This plot shows the CR3BP period of the 
solution vs. the Jacobi integral. For comparison with Figure 26, the highest value of Jacobi integral for the Planar family 
corresponds to the solution with the largest perigee radius value. For the Reflection and Axial families the period varies 
slowly with Jacobi integral. 



Branch Point from 
Planar to Reflection 

Branch Point from 
Planar to Axial 

CR3BP x 

-0.7821 

-0.8628 

CR3BP y 

0 

0 

CR3BP z 

0.0000 

0 

CR3BP v x 

0 

0 

CR3BP v y 

0.0526 

0.3274 

CR3BP v 7 

0 

0.0000 

CR3BP period 

6.2607 

6.2794 

Period (days) 

27.187 

27.268 

Jacobi Constant 

3.1887 

2.9729 

Initial Ra (Re) 

46.40 

51.27 

Initial Rp (deg) 

11.64 

6.85 

Initial inclination (deg) 

0.00 

0.00 

Floquet Multiplier 1 

1.000 

1.000 

Floquet Multiplier 2 

1.000 

1.000 

Floquet Multiplier 3 

1.000 

1.000 

Floquet Multiplier 4 

1.000 

1.000 

Floquet Multiplier 5 

0.7531 +0.6578 i 

0.7023 + 0.7119 i 

Floquet Multiplier 6 

0.7531 -0.6578 i 

0.7023 -0.7119 i 


Table 2. Approximate initial states for the bifurcation points where the Reflection and the Axial families branch from the 
Planar family. The set of Floquet multipliers for each case has two pairs of values +1, showing that each solution belongs 
to two families of periodic orbits. 
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3:1 Periodic Families: Maximum Floquet Multiplier vs. Perigee Radius 



3:1 Periodic Families: Maximum Floquet Multiplier vs. Inclination 



Figure 28. Maximum Floquet multiplier magnitudes for the periodic solutions. The left plot shows the maximum 
multiplier vs. perigee radius. The right plot shows the maximum multiplier vs. inclination for Reflection and Axial 
families. In each plot the red dotted line for the Reflection family was plotted slightly above its actual height of 1 to make 
it more visible. The fact that the maximum Floquet multiplier for the Axial branch approaches 1 for increasing perigee 
radius and inclination indicates that the branch is approaching another bifurcation point. 


We have noted that the Planar Mirror family and the entire Axial family is spectrally unstable. However, 
as shown in Figure 28, even for the spectrally unstable solutions the maximum Floquet multiplier does 
not rise above about 1.017 per lunar cycle for the Planar or Axial families, so in twenty years a 
perturbation would not grow by more than a factor of 1.017 267,6 = 91. Moreover, the orbit where the 
Axial branch becomes vertical appears to be still another bifurcation point, as noted in the caption to 
Figure 28. 

The computation of families of periodic orbits provides a variety of starting points for a mission trajectory 
design. For each periodic orbit we can use properties such as orbit shape, orbit period, the apsidal rotation 
in selecting an orbit. One can also use the linear analysis to identify the quasi-periodic solutions around a 
given periodic solution. 

It is intriguing that the bifurcation diagram for families of nearly 3:1 resonant solutions in the CR3BP has 
some structure in common with bifurcation diagram for some families of libration point orbits shown in 
[14] and [25]. The Reflection resonant orbits exhibit the same mirror symmetry as the libration point Halo 
orbits; the Axial resonant orbits exhibit the same axial symmetry as the libration point Axial or “Y” 
libration point orbits; and the Planar Mirror resonant orbits share both symmetries with the Lyapunov 
libration point orbits. Moreover, the libration point orbits have a bifurcation diagram similar to the one 
shown above for the resonant orbits: The families of Halo orbits and the Axial / “Y” orbits each branch 
from the Lyapunov orbits at different points. The analogy between the two bifurcation diagrams may shed 
further light on both sets of families of periodic solutions. 


The Lidov-Kozai Mechanism 

We observed above that the quasi-periodic oscillations in the IBEX orbit, as well as those around the 
Reflection periodic orbit, exhibit two time scales: one on the order of 10 months, and a second on the 
order of 17 years. It appears that the long-period oscillation is connected with the Kozai resonance 
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mechanism, also called the Lidov-Kozai mechanism [40-41]. The Lidov-Kozai mechanism is believed to 
play a significant role in the evolution of extrasolar planetary systems [26] and triple star systems [27]. 

Ely [28] noted that the Lidov-Kozai mechanism also describes the Earth gravitational perturbation on 
lunar orbiters. The derivation of the Lidov-Kozai model begins by averaging out short-term oscillations. 
The model predicts that both the semimajor axis a and the quantity 

(7) H k = Vl -e 2 cos(i) 

are constant. The quantity yfaH K is one of the Delaunay variables and represents the z component of 
orbital angular momentum [13] [13]. The preservation of H K is also connected with the Tisserand 
criterion in equation (4), since the semimajor axis is constant. It follows that, for prograde orbits, if 
eccentricity increases then inclination decreases, and vice versa. Put another way, the inclination and 
perigee radius vary in unison. We observe this behavior in Figure 17 and Figure 18 for the Reflection 
periodic solution. (Even for the unstable Planar Mirror and Axial solutions, Figure 15 and Figure 21 show 
that the perigee radius and inclination vary together.) Mazeh and Shaman [27] give an approximate period 
of the Lidov-Kozai oscillation: 



In expression (8), Prepresents period, m represents mass and a represents semimajor axis, while for our 
application the subscripts 1, 2 and 3 represent the Earth, spacecraft and Moon, respectively. To derive 
expression (8), terms on the order of unity were neglected, so (8) represents an order-of-magnitude 
approximation. Entering the relevant values in expression (8) we obtain an approximate Kozai period of 
18.5 years, which agrees with the oscillation period for mode Q2 of 16.9 years found for the Reflection 
orbit in Figure 17 and Figure 18. Moreover, for all Reflection orbits with inclinations from 15 deg to 67 
deg, the mode Q2 oscillation period ranges from 37 years down to 9 years. These periods are within a 
factor of 2 of the approximation (8). This strongly suggests that the Lidov-Kozai mechanism is connected 
with the long-term oscillations of the Reflection periodic orbits. 

However, care must be taken when attempting to apply the Lidov-Kozai mechanism to near-resonance 
orbits such as IBEX. The derivation of the Lidov-Kozai model assumes that one can average out short- 
term oscillations. That assumption is not valid near a Mean Motion Resonance, such as the 3:1 resonance 
of IBEX. The Lidov-Kozai mechanism predicts qualitatively different evolution for the eccentricity and 
argument of perigee, depending on whether inclination is above or below a certain critical inclination that 
depends on the semimajor axis. Those predictions are not consistent with our bifurcation analysis, where 
stability of the Reflection and Axial solutions does not change with inclination. Moreover, the Lidov- 
Kozai mechanism does not apply directly to the IBEX orbit because IBEX is significantly affected by the 
Sun as well as the Moon. Thus while Figure 4 and Figure 5 show that the IBEX perigee radius and 
ecliptic inclination generally trend together, the two plots do not vary in unison. There have been studies 
of asteroid motion, such as [29] by Moons and Morbidelli, that consider secular resonances for asteroids 
near the 3:1 resonance in the context of a four-body problem (Sun-Jupiter-Saturn-asteroid). Such analysis 
may shed light on the spacecraft orbital dynamics near mean motion resonance in the Earth-Moon-Sun- 
spacecraft problem. 
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Propagation of an IBEX Orbit State in the CR3BP 


We began this study to gain some understanding of the IBEX orbit dynamics exhibited in Figures 1 to 5. 
We have examined the dynamics of similar, near-resonant orbits in the CR3BP because it can offer 
insights through a simplified system that only models the gravitational effects of the Earth and the Moon, 
although the CR3BP cannot capture the full dynamics. To see how well the CR3BP dynamics capture the 
actual orbit dynamics of IBEX, it is useful to take an initial state from the IBEX orbit, transform it into 
CR3BP coordinates and propagate it using CR3BP dynamics. We can then compare the long-term 
behavior to that seen in Figures 1 to 5 for a high-fidelity force model. At 23 April 2013 00:00:00 UTC, 
the J2000 position of IBEX was (242993.633409, -34985.396684, 52286.448539) km and the velocity 
was (-0.507391, 0.883484, 0.094170) km/sec. Using the state of the Moon at that epoch, we transformed 
this initial state to obtain a CR3BP position (-0.6546, -0.0126, 0.1162) and velocity (0.5619, -0.1561, - 
0.1723). 

Figure 29 shows the evolution of the semimajor axis under CR3BP dynamics. Figure 29 is comparable to 
Figure 3, which shows the evolution in a high-fidelity force model. We see that the semimajor in both 
plots oscillates about a constant value of about 28.9 Earth radii. The range of semimajor axis values in 
Figure 29 is similar to the range in Figure 3, but slightly smaller. The oscillation in Figure 29 exhibits a 
period of about 9 months, close to but slightly longer than the oscillation in Figure 3. In addition, Figure 3 
exhibits more variability that Figure 29, which is to be expected since Figure 3 includes the effects of the 
Sun’s gravity and other forces. 


CR3BP propagation of IBEX semi major axis for 20 years 



Figure 29. Twenty year evolution of IBEX semimajor axis under CR3BP dynamics. The range of values is comparable to those 
in Figure 3, but slightly smaller. The plot exhibits an oscillation with period near 9 months, slightly longer than the oscillation 
seen in Figure 3. 

Figures 30 and 31 show the evolution of the perigee radius and inclination from the lunar orbit plane, 
respectively, under CR3BP dynamics. It is apparent that the perigee radius and inclination oscillate in 
unison, with a period of about 24 years, consistent with the Fidov-Kozai mechanism. Figures 30 and 31 
can be compared to Figures 4 and 5, respectively, for the high-fidelity force model. The CR3BP model 


Page 28 


Perigee Radius (Earth radii 


captures some of the same dynamic features as the high-fidelity model: a medium-term oscillation with 
period around 9 months, and a long-term oscillation with a period of decades. 


CR3BP propagation of IBEX perigee radius for 20 years 



Figure 30. Twenty year evolution of IBEX perigee radius under CR3BP dynamics. 


CR3BP propagation of IBEX inclination for 20 years 
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Figure 31. Twenty year evolution of IBEX inclination from the lunar orbit plane, under CR3BP dynamics. 

However, it is also apparent that the two models show remarkably different long-term evolution. In 
particular, thirteen years after the epoch, the high-fidelity model predicts perigee radius and inclination 
will be near minimum values of about 6 Re and 1 1 deg, respectively, whereas the CR3BP model predicts 
perigee radius and inclination will be near maximum values of about 13 RE and 26 deg, respectively. The 
key factor that is missing from the CR3BP dynamics is the gravity of the Sun. The Sun’s gravity affects 
the evolution of inclination from the lunar orbit plane is two ways: There is a direct effect on the motion 
of the satellite itself, and an indirect effect as the Sun’s gravity causes the lunar orbit plane to evolve with 
a period of 18.6 years. (See [30], p. 281.) In an analysis of the dynamics of the TESS mission orbit near 
2:1 resonance [31], we demonstrated that by including the Sun’s gravity in a Bi-circular restricted 4-body 
problem, we can achieve much better agreement with the high-fidelity propagation. However such a 
model lacks the elegance of the CR3BP that makes it so well suited as a starting point in the design of 
mission where multibody dynamics are important. Time appears explicitly in the dynamics of the 4-body 
problem, and we generally cannot compute continuous families of periodic orbits. 


Conclusions and Future Work 

The goal of this study was to understand better the dynamics of the IBEX extended mission orbit and its 
apparent long-term stability. Toward this end we have examined, in the context of the Circular Restricted 
3-Body Problem, three types of near 3:1 resonant periodic orbits in the Earth-Moon system that are 
similar to the IBEX orbit: Planar Mirror, Reflection and Axial orbits. The results are summarized in Table 
3. 

We first studied these periodic solutions using Floquet theory to assess local stability, and then extended 
the analysis with long-term propagations of perturbations along each of the Floquet modes. The 
Reflection solution is spectrally stable. Although the Planar Mirror solution and Axial solution are 
spectrally unstable, the instability is so mild that it has little practical effect over twenty years. The Lidov- 
Kozai mechanism helps to describe the long-term quasi-periodic oscillations observed in the orbit 
elements, but that model does not consider mean motion resonance and so does not apply completely in 
our study. The Poincare maps showed that all three solutions appear stable for twenty years under 
perturbations more than a thousand times larger than the IBEX orbit determination errors. This makes 
near-resonant orbits very attractive for space missions using small spacecraft with limited propellant 
budgets. 

We also generated bifurcation diagrams of these families of solutions. The bifurcation diagram provides a 
broader understanding of the near-resonant solutions, and offers the mission designer a wider range of 
choices. Using one of the periodic solutions shown above, Randy Paffenroth of Numerica Corporation 
has used the AUTO numerical continuation software to compute a more complete bifurcation diagram for 
the near-resonant families, and has found further bifurcations [32]. 

We concluded the analysis by examining the evolution of an IBEX initial state under CR3BP dynamics. 

In the CR3BP model, semimajor axis in Figure 29 oscillates with period near 9 months and shows no 
secular growth, similar to the behavior for the high-fidelity force model in Figure 3. Perigee radius and 
inclination in Figure 30 and Figure 31 show oscillations with period near 9 months and near 24 years. 
Perigee radius and inclination vary in unison with large amplitude in the long-period oscillation, 
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consistent with the Lidov-Kozai mechanism. Perigee radius and inclination evolution under the high- 
fidelity model in Figure 4 and Figure 5 show more complex features, due to the Sun’s gravity on the 
spacecraft and the Moon. 



Planar Mirror 

Reflection 

Axial 

Analogous Libration Point 
Orbit family 

Lyapunov 

Halo 

Axial 

Design parameter used, 
comparable to IBEX orbit 

Apogee radius = 49.95 Re 

Inclination to lunar orbit 
plane = 23 deg 

Inclination to lunar orbit 
plane = 23 deg 

Floquet theory analysis to 
determine multipliers and 
modes (Figures 6 to 12, and 
Table 1) 

One quasiperiodic pair, with 
period 222 days 

Unstable component grows 
by factor of 10 in 8090 days 

One quasiperiodic pair, with 
period 272 days 

One quasiperiodic pair, with 
period 6167 days 

One quasiperiodic pair, with 
period 222 days 

Unstable component grows 
by factor of 10 in 6719 days 

Floquet mode propagation, 
with perturbation 
magnitude comparable to 
orbit determination error 
(Figures 13 to 22) 

Semimajor axis shows 
quasiperiodic motion but 
very slow growth 

Perigee radius and inclination 
show quasiperiodic and 
stable/unstable components 
with time scales matching 
Floquet multipliers 

Semimajor axis shows 
quasiperiodic motion 

Perigee radius and inclination 
show quasiperiodic with time 
scales matching Floquet 
multipliers 

Semimajor axis shows 
quasiperiodic motion but 
very slow growth 

Perigee radius and inclination 
show quasiperiodic and 
stable/unstable components 
with time scales matching 
Floquet multipliers 

Poincare section analysis, 
with perturbation 
magnitude comparable to 
orbit determination error 
(Figures 23 to 25) 

Twenty year propagation 
with perturbations on order 
of orbit determination error 
demonstrates effective 
stability, and quasiperiodic 
oscillations 

Twenty year propagation 
with perturbations on order 
of orbit determination error 
demonstrates effective 
stability, and quasiperiodic 
oscillations 

Twenty year propagation 
with perturbations on order 
of orbit determination error 
demonstrates effective 
stability, and quasiperiodic 
oscillations 

Bifurcation analysis (Figures 
26 to 28, and Table 2) 

Orbits are mildly unstable for 
perigee radius between 6.85 
and 11.64 Re, otherwise 
spectrally stable 

All orbits are spectrally 
stable. Family branches from 
Planar Mirror family at 
perigee radius 11.64 Re, then 
perigee radius decreases with 
increasing inclination 

All orbits are mildly unstable. 
Family branch from Planar 
Mirror family at perigee 
radius 6.85 Re, then perigee 
radius increases with 
increasing inclination until 
vertical circular orbit in 
inertial frame 

Evolution of IBEX state 
under CR3BP dynamics 
(Figures 29 to 31) 

Semimajor axis shows no secular growth, and oscillates with period of about 9 months, similar 
to high-fidelity dynamics. Perigee radius and inclination show two oscillations with periods 
near 9 months and 24 years. Perigee radius and inclination oscillate in unison, with large 
amplitude, consistent with Lidov-Kozai mechanism. High-fidelity dynamics exhibit more 
complex features, due to Sun's gravity on the spacecraft and on the Moon. 

Practical stability 

Over the course of a decade-long mission, it would be difficult to distinguish long-term 
oscillations from slow-growth instability. Either case would produce practical orbit stability, 
with few or no orbit maintenance maneuvers required. 


Table 3. Summary of analysis results for the three families of resonant orbits in the CR3BP, and for IBEX orbit evolution under 
CR3BP dynamics and high-fidelity dynamics. 


The CR3BP can be a useful tool for understanding the dynamics of a resonant orbit. It allows us to 
capture the qualitative evolution of semimajor axis, perigee radius and inclination, including the 
oscillation of the least two in unison. The CR3BP also allows us to compute families of periodic resonant 
orbits, much as the CR3BP has been used to compute families of Libration Point Orbits. However we 
have seen that the CR3BP does not capture the orbit evolution accurately, primarily because it neglects 
the Sun’s gravity on the spacecraft and on the Moon. The TESS mission orbit will be near 2: 1 resonance 
with the Moon. In a study of the TESS orbit similar to this one [31], we included the Sun’s gravity using a 
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Bi-circular Restricted 4-Body Problem (BCR4BP), and achieved good agreement with the high-fidelity 
force model for the duration of a four-year mission. 

There are several directions for future work in this area. Like our study of the TESS orbit, we can use the 
BCR4BP to study other resonant orbits. The present study has focused on periodic solutions, but it is clear 
that the IBEX orbit is a quasiperiodic solution, so quasi-periodicity deserves further attention. We have 
not yet addressed how to reach a resonant orbit, but Vaquero and Howell have investigated transfer 
between resonance orbits via invariant manifolds [27-29]. Works such as [29] on asteroid motion may 
help us to understand better spacecraft orbit dynamics near a mean motion resonance. In particular, while 
the Lidov-Kozai mechanism appears to describe the long-term oscillations in these orbits, we do not yet 
have a model for the oscillations with period near 9 months. 

The present study has focused on 3:1 resonant solutions. It would be interesting to examine orbits near 
2:1 and 4:1 resonance, both of which we considered for the IBEX extended mission orbit [1] and may be 
useful in future missions. In [31] we have generated a bifurcation diagram for the 2:1 periodic orbits, and 
found that near the Planar family the bifurcation diagram for the 2:1 case is very similar to the 3:1 case. 
However the global behavior of the 2:1 Axial branch appears distinctly different from the 3:1 case: We 
did not find a bifurcation point in the 2:1 Axial family until the inclination reached 180 deg, where the 
orbit is circular. Finally, it would be intriguing to pursue the analogy between resonance orbits and 
libration points orbits in the CR3BP. 
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Appendix 


The CR3BP system (1) can be written concisely as 

(dH\ T 

(ad s = m = K (-) 

T 

qh T fdH dH\ 

where — = is the derivative of the Hamiltonian H(s) = —C(s)/ 2, and C(s ) is the Jacobi 

integral given by equation (2). The 6-by-6 skew-symmetric matrix K is given by 



where 0 is the 3-by-3 zero matrix, I is the 3-by-3 identity matrix, and 


h = 



1 

0 

0 



Equation (A.l) is a nonstandard Hamiltonian system. In a standard (or canonical) Hamiltonian system, the 
block matrix in the lower right of matrix K would be 0. It is also possible to express the CR3BP as a 
standard Hamiltonian system, where the velocity is replaced by momentum variables in the state vector 
[ 11 ]. 


Let s(t, 5 0 ) be the solution of equation (A.l) with state s 0 at time t — 0. Let 0(t, s 0 ) be the state 
transition matrix (or fundamental matrix) for system (A.l), so that 0(t, s 0 ) is the derivative of the 
solution with respect to the initial state. That is, (t, s 0 ) = ds(t, z)/dz | z=5o . The matrix <J>(£, s 0 ) 
satisfies the differential equation 

(A. 2 ) <D' = Df(s(t,s 0 ))<I>, <KO,s 0 )=I 
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Equations (1) and (A.2) together define a system of 42 first-order differential equations for state s and 
state transition matrix O. We sometimes suppress the dependence of O on the initial state for brevity. 
Because the system (A.l) is Hamiltonian, 

(A3) Df(s) = KD 2 H(s) 

where D 2 H(s)= d 2 1 1 (s) / d 2 S is a symmetric matrix. 

The state transition matrix satisfies the symplectic condition 
(A. 4) 0 t K“ 1 0= K 1 

Equation (A.4) is true at t = 0 because (0, s 0 ) = I . To complete the proof we can differentiate 
<j> t K - 1 <J> with respect to time, then employ equations (A.2) and (A.3) and the fact that K is skew- 
symmetric. One implication of equation (A.4) is that det(0(t)) = +1 for all t, which implies that the 
product of eigenvalues of O equals +1. Because 4> is invertible, 0 is not an eigenvalue. Suppose 
<J> v = 2 v where v =£ 0. Because O is real-valued, O v = X v so the conjugate X is also an eigenvalue. 
If we pre-multiply O v = A v by O t K _ 1 , it follows that (K -1 v) T is a left eigenvector of O with 
eigenvalue 1/2. It follows that eigenvalues of O other than ±1 occur in groups: 

• A pair (A, 1/2) where 2 is real and 2 =£ 1 

• A pair (A, 2) where 2 is complex and |2| = 1 

• A quadruplet (A, A, 1/2 1/2) where 2 is complex and |2| ^ 1. 

The first two of these cases are the most common. 

Let s p (t, Sq) be a periodic solution of (1) with period T , and let M = 0(T) be the corresponding 
monodromy matrix. From Floquet theory, the state transition matrix can be expressed as 

(A. 5) <D(t - t 0 ) = P(t)exp(B(t - t^PCt,,)- 1 

where P(t) is nonsingular with period T , and B is a constant matrix. See, for example, [1], or [18]. The 
monodromy matrix M is then given by M = exp(B T ). Equation (A.5) shows that the evolution of a 
perturbation is governed by two distinct factors: The constant matrix B in the exponent describes the 
long-term variation, and the periodic matrix P(t) describes the oscillation within a single period. For 
stability analysis, we are only concerned with the long-term behavior. Floquet theory implies that, to 
assess the linear stability of a periodic solution, it is sufficient to examine the eigenstructure of the 
monodromy matrix M. As shown in [11] (Lemma 6.4.1, p. 77), because the dynamics (1) are time- 
independent, the state derivative /(Sq) i s a right eigenvector of M with eigenvalue +1: 

(A 6) M f(sS)=f(sS) 

In addition, dH/ds is a left eigenvector of M with eigenvalue +1. This fact is proved in [1 1] (Lemma 
6.5.1, p. 80.). It also follows because f(so) = K (5H/ ds) T is a right eigenvector with eigenvalue +1, so 

^K _1 f(so)i = dH/ 9s is a left eigenvector with eigenvalue +1. 
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In a family of periodic solutions, the period generally depends on an integral of motion. In the case of 
Keplerian motion, the period of an orbit depends on the energy, but not on the angular momentum or the 
Runge-Lenz vector. Figure 27 shows that period depends on the Jacobi integral in each of the three 
families studied here. Wiesel [18], p. 84, notes that, if the period depends on an integral C, then M must 
have a generalized eigenvector associated with eigenvalue +1: If v be a vector along /(Sq), then there 
exists a generalized eigenvector g and a scalar r such that 

(A 7) M (/5) = (/5) (J l) 
or 

(A8) M / = / 

( 71 . 9 ) Mg = g + rf 

Moreover, a perturbation along g perturbs the integral C and so perturbs the period. Equation (A. 8) is a 
restatement of equation (A.6). Equation (A.9) also has a basic physical interpretation. Suppose that the 
initial state of the nominal periodic solution Sq is perturbed a small amount along direction g so that the 
integral C is perturbed slightly from the nominal. Then the period of the perturbed solution is also slightly 
different from the nominal solution, so after one period the perturbed solution will be offset by a small 

amount along the tangent direction. As noted above, this is precisely the behavior we see when a 

Keplerian orbit is perturbed in a way that changes the energy. 

When we numerically compute the eigenvalue decomposition of monodromy matrix M using the Matlab 
function eig , the eigenvalue +1 cannot be found precisely, due to finite precision. Instead Matlab can 
produce a pair of eigenvalues 1 + a and (1 + £) _1 « 1 — £, with corresponding eigenvectors of the form 
/ + y 9 and f — y g, where £ « 1 and y « 1. Therefore 

(A. 10) M(/ + Y0) = (1 + £ )(/ + 75)= / + £ f+ 7 5 + £ Y9 

(A. 11) M (f - y g) = (1 - e)(f - 7 5)= f~ £ f~ 75 + £ Y9 

If we neglect terms of order ey in equations (A. 10) and (A.l 1), and take the sum and the difference of 
these two equations, we obtain equations of the form (A.8) and (A.9), where = e/y . Thus we can use the 
Matlab results to construct the generalized eigenvector. Rather than producing a pair of real eigenvalues 
near 1, Matlab eig can produce a complex conjugate pair of eigenvalues near +1. In that case, we can 
separate the real and imaginary pairs of the eigenvectors to obtain vectors / and g. 
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